Optimal ramp shapes for the fermionic Hubbard model in infinite dimensions 



Nikolai Eurich, 1 Martin Eckstein, 1 and Philipp Werner 1 

' ' Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland 
(Dated: May 3, 2011) 

We use non-equilibrium dynamical mean field theory and a real-time diagrammatic impurity solver to study 
the heating associated with time-dependent changes of the interaction in a fermionic Hubbard model. Optimal 
ramp shapes U (t) which minimize the excitation energy are determined for a noninteracting initial state and an 
infinitesimal change of the interaction strength. For ramp times of a few inverse hoppings, these optimal U(t) 
are strongly oscillating with a frequency determined by the bandwidth. We show that the scaled versions of 
the optimized ramps yield substantially lower temperatures than linear ramps even for final interaction values 
comparable to the bandwidth. The relaxation of the system after the ramp and its dependence on the ramp shape 
are also addressed. 
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I. INTRODUCTION 

Correlated electron systems play a central role in mod- 
ern condensed matter physics. 1 While most theoretical studies 
have focussed on the equilibrium properties of these systems 
and many interesting phenomena remain to be understood, 
there is a growing effort to also address the non-equilibrium 
dynamics. These investigations are motivated by ongoing ex- 
perimental progress in this field, in particular measurements 
of the relaxation dynamics in correlated materials by means 
of pump-probe spectroscopy, 2 and experiments on cold-atom 
systems^ The latter provide an experimental realization of 
theoretical models such as the fermionic one-band Hubbard 
model,~ and allow to explore the time-evolution after rapid 
changes in model parameters^ 

A computationally tractable method to simulate the time 
evolution of a correlated lattice model after a quench or ex- 
ternal perturbation is nonequilibrium dynamical mean field 
theory, 7,8 This approach neglects spatial correlations, and lo- 
cal correlation functions and their dynamics are obtained from 
an appropriately defined quantum impurity model. The re- 
peated numerical solution of this impurity model on the real- 
time (Keldysh) contour is the most challenging aspect of 
nonequilibrium DMFT calculations, but several perturbative 9 
and exact 10 impurity solvers have recently been developed to 
tackle this problem. Nonequilibrium dynamical mean field 
theory has been applied to interaction quenches in the Hub- 
bard model 12-14 as well as in the Falicov-Kimball model^ 
and to study the time evolution in the presence of external 
fieldsi^r— Interesting phenomena have been observed in 
these theoretical investigations, such as an apparent dynamical 
phase transition in the relaxation dynamics after an interaction 
quenchJ 2 *^ 

A change in the interaction parameter affects the total en- 
ergy and therefore the temperature of the system after equili- 
bration, and nonadiabatic heating effects can be substantial. 20 
The dynamical phase transition after a sudden interaction 
quench in the T = 0, noninteracting Hubbard models 2 , for 
example, occurs at an energy which translates into a temper- 
ature of about 0.2 times the bandwidth, which is far higher 
than the temperatures for which a metal-insulator transition is 
found in equilibrium. While it is an interesting observation 



that an apparently sharp dynamical transition occurs in such 
a highly excited system, one would like to avoid the heating 
effect in other contexts. In cold-atom systems, where much 
effort is devoted to realize a fermionic Hubbard model in the 
low-temperature regime, interactions and the optical lattice 
must be switched on in such a way that the heating is min- 
imized. If there are constraints on the ramping time, it can 
become a subtle question whether an equilibrium state in a 
desired low-temperature phase can be reached at all by means 
of a ramp which starts from another, more easily preparable 
state of the system. 21 And since in most cases one will be in- 
terested in preparing an equilibrium state of the interacting 
system, a second relevant question concerns the time-scale 
on which the system relaxes and its dependence on the ramp 
shape. Another example is the recently proposed conversion 
of repulsive interactions into attractive ones by means of ex- 
ternal periodic electric fields (also essentially an interaction 
quench), 19 which raises the interesting possibility of AC-field 
induced superconductivity. Also in this context the heating 
and the thermalization time associated with the ramp protocol 
are important issues. 

The purpose of this theoretical study is to determine the 
optimal ramping procedure between two different parameter 
regimes of the Hubbard model, where "optimal" means the 
passage in which the system is least excitedi 2 ^ 2 ^ In particular 
when the ramping time is restricted to only a few times the 
inverse hopping, one can expect a strong dependence of the 
excitation energy and of the relaxation dynamics on the ramp- 
ing protocol. However, an unbiased optimization over the in- 
finite space of ramp shapes is not possible, because even the 
computation of the excitation energy for a single ramp proto- 
col requires considerable numerical effort. For "small ramps", 
on the other hand, in which the parameter is changed by only 
a small amount, one can resort to perturbation theory in the 
ramp amplitude to disentangle the influence of the ramp shape 
and many-body effects on the excitation energy^ 4 This allows 
us to perform an efficient optimization over a wider space of 
ramp shapes. In the present work we demonstrate that ramp 
shapes which are optimized for such infinitesimal parameter 
changes yield a considerably lower excitation energy than the 
generic linear ramp or other simple ramping protocols, even 
when applied to ramps with an amplitude comparable to the 
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FIG. 1: Interaction ramp from U(t = 0) = to U(t 
simulated time is t max > r. 
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bandwidth. We will furthermore demonstrate that a dynamical 
transition (associated with fast thermalization) also exists in 
the case of short ramps, similar to what has been found for an 
interaction quench This observation suggests that a suit- 
able choice of the ramp shape allows the preparation of ther- 
mal equilibrium states over a wide range of interaction values, 
within a switching time of only a few inverse hoppings. 

The specific model we consider is the one-band Hubbard 
model, 



2> (1) 



with hopping amplitudes and a time-dependent on-site re- 
pulsion U(t). The hoppings Vij are chosen corresponding to 
a semi-elliptic density of states of bandwidth 4V, p(e) — 
V ' W 1 — € 2 /(2ttV 2 ) and we restrict our calculations to the 
paramagnetic phase of the half-filled model. We set V = 1 
as the unit of energy. Initially, the system is prepared in the 
noninteracting ground state (U(t = 0) = 0, T(t = 0) = 0). 
The interaction U(i) is switched from to U in a time r and 
then kept fixed at U for t > t (see illustration in Fig. [TJ. We 
explore different ramp shapes for switching on the interaction 
U(t) and compute the resulting heating effect, which is quan- 
tified by the excitation energy. Since the Hamiltonian is time- 
independent for t > r, the excitation energy is constant for 
t > r, even though other observables take time to equilibrate 
after the end of the ramp. 

The rest of the paper is organized as follows. In Section III1 
we briefly describe the real-time Monte Carlo impurity solver 
and the small adaptations needed to treat time-dependent in- 
teractions, in Section UTU we recall the results of Ref. for 
the infinitesimal quench and show how U (t) can be optimized 
based on these second order perturbation theory formulas. 
In Section [TV] we present nonequilibrium DMFT results for 
U = 1 and U = 3 and compare the optimized shapes to linear 
ramps. We will also discuss the time evolution of the momen- 
tum distribution function after the ramp and demonstrate that 
a fast thermalization occurs at a well-defined, but ramp shape 
dependent value of the interaction. Section [V] is an outlook 
and conclusion. 



II. WEAK-COUPLING QUANTUM MONTE CARLO 

We compute the real-time evolution of the impurity model 
within the DMFT selfconsistent loop using a continuous-time 
Quantum Monte Carlo (CT-QMC) method^ which is based 
on an expansion of the partition function in powers of the 
interaction term. 25 This method is free of systematic errors 
but restricted to relatively short times due to a dynamical 
sign problem. Nevertheless, it has proven useful for interac- 
tion quench calculation a 12 ' 13 ' 19 in the Hubbard model and for 
the simulation of transport through quantum dots>2& Here, we 
briefly recapitulate the weak-coupling formalism to show that 
it can account for time-dependent interactions, simply by re- 
placing weight factors associated with individual vertices by 
time-dependent factors. The remaining part of the DMFT self- 
consistency, i.e., the computation of lattice observables and 
lattice Green functions from the local Green function of the 
impurity system, is unchanged with respect to the interaction 
quench setup. A detailed description of these equations can 
be found in Ref. [TH 

We start by writing the partition function Z = Tr[e~P H ] 

as 



Z = Tr 



-P Hi 



Tee 



- f c dtH 2 



(2) 



where J c denotes the integral along the Keldysh-contour — > 
tmax — > — > —ij3. In the simulations, we choose the length of 
the contour t max somewhat larger than the ramp time r. The 
operators Hi and H2 correspond to the hopping and interac- 
tion part of the impurity Hamiltonian. Specifically, we use 



Hi 



H v - k(t)/U 



with 



H u = U(t){n t n ir -(n t + n l )/2). 



(3) 



(4) 



The nonzero real function k(t) was introduced to enable an 
auxiliary field decomposition of the interaction term. Expan- 
sion of the contour-ordered exponential in powers of Hjj — 
k(t) /t max and the application of the decoupling formula^ 



l-t max Hu/k(t) = - ]T 

S- 

cosh(7(t)) = 1 + 



= 7(t)s(n t - nj.) 



s=-l,l 

2k(t) 



(5) 



(6) 



at every interaction vertex results in an expression of the par- 
tition function as a sum over all possible collections of Ising 
spin configurations on the contour with weight 

u({(ti,si), (t 2 ,s 2 ), (t n ,s n )}) = 

(-i n -)(i n +)(k(t)dt/(2t max )) n - +n +l[detN- 1 . (7) 

a 

Here, n + and n_ are the number of spins on the forward and 
backward branch of the contour, and we have used the fact 
that in a noninteracting initial state there are no spins on the 
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imaginary-time branch. The matrices N^ 1 are determined by 
the location of the interaction vertices on the contour C and 
are given by 



N- 1 = e s ° 



-(iG ,a)(e s ° -I), 



(8) 



with Go.cr the bath Green's function and e f 



diag(e 



7(ti)si<T 



., e 



In the actual calculations, we 



choose j(t) = 7 time-independent, so that the time depen- 
dence of the interaction manifests itself only in a time depen- 
dence of the parameter k(t) = ^t miix U (t) / (cosh^) — 1). 

The sampling procedure consists of generating all spin con- 
figurations on the contour C through random insertions and 
removals of spins. During the sampling, we measure the quan- 
tity 

X a (si,s 2 ) = 

n 

V M*i^)[(e 5 " - lW^Scfatj)) .(9) 
* — ' ' MC 



The impurity Green's function is then obtained as 

G a (ti,t 2 ) = Go t<7 (ti,t2) 

+ J dsi J ds 2 G , a (ti,si)X a (si, s 2 )G ^{s 2 ,t 2 ). (10) 
m. PERTUBATIVE ANALYSIS FOR SMALL RAMPS 
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A. General considerations 



In this section we briefly recall the perturbative analysis of 
Ref. [24] for an infinitesimal interaction ramp. To be specific, 
we consider a ramp of the form 



U(t) = Ur(t/T), 



(ID 



where r(x) is the ramp-shape function which satisfies r(x) = 
for x < and r(x) = 1 for x > 1. The excitation energy 
per lattice site is defined as 



AE(r) - E(t) - Eo(t), 



(12) 



where Eq (r) is the groundstate energy per lattice site of the 
model with interaction U, and E(t) = (H(t)) is the energy 
expectation value of the system with time-dependent interac- 
tion. After expanding in powers of U, the second-order result 
for the excitation energy can be decomposed into contribu- 
tions from a ramp spectrum F(x), which depends only on the 
ramp shape but not on the properties of the system, and from 
an excitation density R(u>), which depends on the system but 
not on the ramp shape,— 



AE(t) = 
£(t) = 

F(x) = 



U 2 £{t) + 0(U 3 ), 

f°° du „, . „. x 
/ —R(u>)F(wt), 
Jo w 

ds r'(s)e v ' 



(13) 
(14) 

(15) 



FIG. 2: Spectral density R{u)/uj and ramp spectra F(ojt) for linear 
ramps with r = 1, 1.25, 2.5 (top panel), as well as ramps of the form 
r(x) = x + 0.2 sin(27rx) and r = 1, 1.25, 2.5 (bottom panel). 

In this expression, the excitation density is defined by the 
Lehmann representation, 

R{u) = ^Y^^ W \^ 5 ^- En + E °^ (16) 

n#0 

where W = (n-j- — \){n^ — |) is the operator which couples 
to the time-dependent parameter U(t) in the Hamiltonian, L 
is number of lattice sites, and \<f> n ) and E n are the eigenfunc- 
tions and eigenvalues of the Hamiltonian. The function R(uS) 
can be evaluated easily for the Hubbard model with [7 = 0, 
leading to^i 



dep{e) 



dnp{fi 



;) / dvp{v) P {[l 

Jo 



-v). 
(17) 



In Fig.|2]we plot the function R(cj) /w for the semi-elliptical 
density of states together with ramp spectra for different ramp 
shapes. At U = 0, the only relevant energy scale is the hop- 
ping V (= 1): R(ui)/uj has a peak near ui ps 3 and vanishes 
for lu > 8 = 2 • bandwidth. According to Eq. ( [Pil l, a small 
overlap of the excitation density and ramp spectrum results 
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in a small excitation energy. If the ramp spectrum F(x) falls 
off rapidly at large x, the main weight of F{ujt) and the main 
contribution to the integral in Eq. (TBI i comes from frequencies 
u < 1/t. Using the asymptotic form R(lj) oc lj 3 for small uj 
one can prove that the system approaches the adiabatic limit 
t — y oo with a power law behavior AE(t) oc 1 /t 3 , provided 
that the high-frequency tail of F(x) falls off faster than 1 /x 3 . 
For the linear ramp r(x) ~ x, for which F(x) falls off as 
1/x 2 , one finds AE(t) cx 1/t 2 ^. 

In the present paper we are interested in the excitation en- 
ergies for ramp times which are so short that the asymptotic 
power law behavior AE(t) oc 1/t 71 does not yet hold, and 
we attempt to minimize AE(t) in Eq. ( fT4l i with respect to 
the ramp shape r{x). The reduction of the excitation en- 
ergy can be understood as a consequence of a suppression 
of the ramp spectrum F(ujt) in the frequency range around 
the maximum of R(oj)/uj, which is achieved by suitably de- 
signing the ramp r(x). Clearly, the optimal ramp shape will 
strongly depend on the ramp time r. For example, for the 
ramp r(x) = x + 0.2 sin(27r.T) the function F(lut) is sup- 
pressed around cj = 3 for r « 1.25, while for t = 2.5 the 
same ramp shape is apparently not favorable (Fig. [2})). 



B. Minimization in a tent basis 

To find the optimal ramp-protocol r{x) which minimizes 
the excitation energy 8(t) given some constraints on r(x), we 
discretize the interval [0, 1] into an equidistant mesh of N + 2 
points x k = kAx(k = 0,...,N+l;Ax = l/(N+l)), and 
consider ramp functions which lineraly interpolate between 
the values r(x k ) = x k + c k . The numbers c k measure the 
deviations from the linear ramp r(x) = x, so Co = c^+i = 0. 
Equivalently, this means that the function r(x) is expanded in 
a basis of tent-shaped functions 4> k (x), 



N 



r(x) = x + y^ c k (j>k(x), 



(18) 



with 



k=l 

\fxe [x k -i,x k ], 

ifxe [x k ,x k+ i], 



■Ck+l- 



Xk + l—Xk 

otherwise 



(19) 



The first derivavtive of the ramp function r(x) is given by 



N 



r'(x) = l + J2 c Mx) 



(20) 



k=l 



with 



3j if x G [xk-i,x k ], 
4>' k {x) = { ^| if x e [x k ,x k+ i], 
otherwise. 



(21) 



The calculation of F(ljt) requires the evaluation of the inte- 
gral 

1 



r'(s)e iujTS ds = - — (e tuJT - l) 



2(1 - cos(wrAx)) 
iujtAx 



N 



(22) 



k=l 



Multiplying the right hand side of Eq. d22l with its complex 
conjugate and remembering that F(x) is real gives 



£(t)= fdu^FH 
= const + c T f + c T Mc, 



(23) 



where the coefficients c k are written as an TV-component vec- 
tor c, and the iV-component vector f and the N x N matrix 
M are given by 



fk = duj 
Jo 



R(u) 4(1 - cos(wrAx)) 



w \ut\ Ax 

x [(cos(cjt(1 - x k )) - cos(u!Tx k )} , (24) 
R(cj) 4(1 - cos(o;TAa;)) 2 



duj- 



kk' i ,9 

o w \u>t\ Ax 2 

x cos(o;r(.T fe -x k >)). 



(25) 



The minimization of this quadratic problem can be performed 
using standard techniques. 



C. Unconstrained optimization and symmetries 

In order to get a qualitative understanding of the optimal 
ramp shapes we first note the following symmetry properties 
of M and f : 

1. M is symmetric, i.e. Mi j = Mj^, 

2. Mi j is constant on each subdiagonal \i — j\ = k, 

3. f is antisymmetric, i.e. /; = — /jv+i-j. 

The optimal, unconstrained ramp shape can be found by 
solving d£(T]{ci})/ dci = and therefore, using the sym- 
metry of M 



2Mc + f = 0. 



(26) 



It follows that the solution c is antisymmetric, i.e. Cj = 
— c/v+i-i- In fact, setting N = 2k, the first i = 1, . . . , k 
components of Eq. d26l i are 

A; 

- fi=/ J ( M i,j c 3 + M it 2k+l-jC2k+l-j) ■ (27) 
J'=l 

Using the symmetry f t = -/ 2 fc+i-i we get 
fc 

= ( M i,j + Mak+i-ij) Cj 



3=1 



(Mi^k+i-j + AI-2k+i-i,2k+i-j) c-2k+i~j- (28) 
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The symmetry properties of M thus imply 

k 

= ^2 i M i,j + Mi,2k+l-j) {Cj + C 2k +l-j), (29) 
i=i 

which can only be satisfied for all i = 1, . . . , N if Cj = 

— C2k-j + l- 

IV. RESULTS 
A. Small ramp amplitude 

We first consider unconstrained paths and solve Eq. ( |26] i 
for different values of r. As illustrated in Fig. [3] the opti- 
mal paths oscillate around the linear ramp, and the number 
of oscillations increases with increasing r. Thus, the some- 
what unexpected result is that paths which minimize the ex- 
citation energy, at least according to the perturbative formula 
(fl4] i. may involve excursions to positive and negative interac- 
tion values which are much larger in absolute value than the 
final interaction U. 

Note that the uncontrained optimization as decribed here 
becomes numerically unstable for large N. This is because 
one can always add to the ramp a highly oscillating component 
Sr(x) whose ramp spectrum 6F(x) lies almost completely 
outside the support of the excitation density R(uj), and thus 
does not influence the excitation energy. In other words, there 
are directions in the parameter space along which the excita- 
tion energy hardly changes. The corresponding eigenvalues 
of the matrix M are almost zero, and the linear equation d26l i 
becomes ill conditioned for large N. 

However, optimal ramps with large amplitude oscillations 
and sign changes may be difficult to realize in experiments. 
We therefore also compute optimal ramps with the constraint 
< r(x) < 1, using the reflective Newton method^ im- 
plemented in MATLAB4i The resulting paths for different r 
are shown in Fig. [4] It is evident from these results that the 
optimal ramp shapes are characterized by a roughly constant 
oscillation frequency u>q = 27rn osc /r w 9 — 10. To gain some 
insight into the origin of these oscillations we plot in Fig. [5J 
the ramp spectra F(uit) for the ramps shown in Fig. [4] The 
inset shows a close-up view of the optimal ramp spectra in the 
frequency range where R(u>) is large. We see that F(uit) is 
optimized in such a way that the overlap with R(u>) is min- 
imal. The main panel shows the ramp spectra in the range 
< u> < 14. A large peak in F(ljt) is evident at cj rj 9.5 
(r = 1.25), 11 (r = 2.25), 9.75 (r = 3.25) and 9 (r = 4.25), 
i.e., just above the largest value of us for which R(ui) > 0. 
We can understand from Eq. ( [TBI that a sharp peak in F(ujt) 
at lo ~ cjo corresponds to an oscillating r'(s) ~ cos(luqts). 
So, the frequency which appears in the optimal ramp shapes is 
determined by the support of the function R(oj), which itself 
is defined by the density of states p(u>). In the case of a sym- 
metric p(uj) considered here, the support of R{oj) is twice the 
bandwidth (Eq. (fTTI)). The optimal ramp shapes have there- 
fore an oscillating component with an oscillation frequency 
roughly given by ojq ~ 2 ■ bandwidth. 
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FIG. 3: Optimized unconstrained ramp shapes for r = 1.25, 2.25, 
3.25 and indicated number of basis functions (N). 
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FIG. 4: Optimized constrained ramp shapes for r = 1.25, 2.25, 
3.25, 4.25 and N = 20. 



B. Monte Carlo results for larger U 

In this subsection, we use Monte Carlo simulations to com- 
pute the heating effect for ramps with amplitudes beyond the 
perturbative regime, and compare the heating produced by lin- 
ear ramps to the heating produced by ramps which were opti- 
mized for infinitesimal ramp amplitudes. 

First, we would like to confirm the validity of the pertur- 
bative analysis for ramps to small interaction. In Fig. [6] we 
plot the excitation energy computed by means of nonequi- 
librium DMFT for different ramps to U = 1. The red line 
with crosses, and the green line with stars show the result 
of Eq. ( TPfl i for linear ramps and ramps of the form r(x) = 
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FIG. 5: Ramp spectra F(lut) for the optimized constrained ramp 
shapes in Fig.g](r = 1.25, 2.25, 3.25, 4.25 and N = 20). The inset 
shows a close-up view of the energy range 1 < w < 6.5 in which 
the excitation density R(cj) is large and hence the ramp spectra are 
strongly suppressed. 



FIG. 7: Effective temperature for linear ramps (black squares) and 
optimized constrained ramps (blue circles) to U — 3 as a function 
of ramp-up time r. The inset shows the relationship between the 
energy after the ramp and the temperature of an equilibrium model 
with interaction U and the same total energy. 



, Eq.(14) 

r(x)=x+0.2sin(27D<), Eq.(14) — >K- 
r(x)=x, QMC o 
r(x)=x+0.2sin(27[x), QMC □ 




FIG. 6: Comparison of the DMFT data (QMC) to the results ob- 
tained from the perturbative analysis (Eq. I14t for both linear and 
sinusoidal ramps to U = 1 and for various ramp times r. The in- 
teracting ground state was estimated to be Eq = 0.6195. 



x + 0.2 sin(27rx), respectively. Blue circles and pink squares 
show the excitation energy obtained from the DMFT calcu- 
lation for an interacting ground state energy Eq = 0.6195. 
This value of Eq gives the best agreement between analytical 
and DMFT results, and is consistent within error bars with the 
ground state energy E^ MFT (U = 1) = 0.618(2) estimated 
from equilibrium DMFT calculations. Hence, for ramps to 
[7 = 1 the formula ( fl~4b gives accurate excitation energies. 
The ramp spectra F(ujt) for several values of r are plotted in 
Fig. |2] and explain the nonmonotonic behavior of the excita- 
tion energy in the case of r(x) = x + 0.2 sin(27rx). 

A nontrivial question is whether the ramp shapes opti- 
mized using Eq. ([Pit yield low excitation energies also for 
larger values of U. In the following we will consider ramps 
to U = 3. This value is close to the critical interaction 



strength U c w 3.2 for instantaneous quenches, where the sys- 
tem was found to thermalize within a time of less than two 
inverse hoppings ^all To quantify the heating effect we trans- 
late the excitation energy into an effective temperature (inset 
of Fig. |7J, which is defined as the temperature for which an 
equilibrium system with interaction U has total energy E(t). 

The effective temperature for linear ramps to U — 3 is plot- 
ted as a function of ramp time in the main panel of Fig. [7] 
(black line with squares). We see that the temperature after 
the quench drops rapidly with r until about r w 1.5 and then 
more slowly for longer ramp times. The times t which are 
accessible with Monte Carlo are not sufficient to determine 
the asymptotic behavior of the excitation energy (r — > oo), 
but it is expected to be AE(t) ~ 1/r 2 based on the pertur- 
bative anaysis. Figure [8] shows the effective temperatures for 
t = 1.25 obtained from the optimized unconstrained ramps 
for different numbers N of basis functions (see top panels of 
Fig. 0. The lowest temperature for r = 1.25, N = 5 and 
t = 2.25, N = 7, are indicated by the red and green crosses in 
Fig. |7] For short ramp times (r = 1.25), the effective tempera- 
ture can be reduced by about a factor of 2 (compared to linear 
ramps) if an oscillating ramp shape is used. For r = 2.25, 
the reduction is only about 10%. The constrained optimized 
ramps yield comparable reductions in the effective tempera- 
ture. In Fig. [7] the blue line with circles shows the effective 
temperature for N = 25 basis functions and the constraint 
< r(x) < 1. 



C. Accuracy of the perturbative results 

The results shown in Fig. UJ demonstrate that the pertur- 
bative analysis of Section [TIT] may be used to compute ramp 
shapes with considerably lower excitation energy than linear 
ramps. On the other hand, it is not yet clear how quantita- 
tively accurate the estimated energies are in the case of ramps 
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FIG. 8: Effective temperatures for optimized unconstrained ramps 
(see top panels of Fig. [3} with N basis functions. U = 3, r = 1.25. 



to intermediate or strong U. If the difference between the 
true and the estimated excitation energy is large and the error 
depends sensitively on the ramp shape, then the true optimal 
ramp might look very different from the shape obtained by our 
procedure. To test the reliability of our estimated excitation 
energy and its dependence on the ramp shape we consider the 
simple one-parameter families r±(x) = x + asin(27ra), and 
7*2 (x) — x + a sin(47rx) with r = 1.25 and [7 = 3, and com- 
pare the optimal values for the parameter a obtained from the 
perturbative analysis to the Monte Carlo results. The oscilla- 
tion frequency of t\ (x) is not compatible with the ramp time, 
while the oscillation frequency of 7-2(2) is identical to that 
found in the constrained optimized infinitesimal ramp (top left 
panel of Fig. [4j and thus results in lower excitation energies. 

In Fig. [9] the blue line with squares plots the excitation en- 
ergy obtained from Eq. ( fT4t as a function of the shape param- 
eter a. The red line with circles shows QMC results for the 
total energy which are shifted in such a way that the minima 
of the curves coincide. The necessary shifts Eq = —0.279 
(for J"i(x)), and Eq = —0.277 (for r^x)) are still compa- 
rable to the ground state energy Eq MFT = -0.284(3) ob- 
tained from aT^O extrapolation of equilibrium DMFT en- 
ergies for the [7 = 3 model, but the perturbative formula ap- 
pears to underestimate the excitation energy of the system by 
8E Q ps 0.005 - 0.007. Still, the variation of the excitation en- 
ergy as a function of the ramp-parameter a is quite accurately 
reproduced. In particular, the perturbative analysis yields the 
correct values for the optimal parameters (a opt = 0.175 for 
r\(x) and a opt = 0.87 for ^(x)), and thus the correct optimal 
ramp shape. 

For interactions U > 3.5, the analytical prediction for the 
excitation energy is no longer in quantitative agreement with 
the Monte Carlo results. In the lower panel of Fig. |9]we show 
Monte Carlo data for U = 4 and 5, with an arbitrary off- 
set (chosen in such a way that the minima are at around the 
same energy). We see that as the interaction is increased, the 
minimum in the excitation energy shifts to smaller values of 
a, while the curvature at the minimum increases. This implies 
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FIG. 9: Excitation energies for the one-parameter model r(x) = 
x + asin(27ri) (top panel), and r(x) — x + asin(47r:r) (bottom 
panel) for r = 1.25, U — 3. The blue line with squares shows 
the result from second order perturbation theory, while the red line 
with circles shows the Monte Carlo data shifted in such a way that 
the minima of the curves coincide. In the lower panel we also plots 
Monte Carlo results (with arbitrary energy offset) for U = 4 and 
(7 = 5. 



that the minimization of Eq. ( fT4b yields ramp shapes which do 
not minimize the excitation energy, and which in the large-C7 
limit may be qualitatively wrong. 



D. Relaxation after the ramp 

We would finally like to address the relaxation dynamics 
after the ramp. If the goal of an experiment is to prepare an 
equilibrium state at given interaction strength, then it is not 
only the ramp time, but also the thermalization time, which 
determines the relevant time scale of this process. The re- 
laxation dynamics after an (sudden) interaction quench has 
been studied in Refs. [l2l - [l4l These calculations demonstrated 
the existence of a "dynamical phase transition" at some crit- 



ical interaction strength U[ 



quench 



3.2, which separates two 



qualitatively different relaxation regimes. After a quench to 



U, 



quench 



, the system thermalizes within a few inverse hoppings, 



whereas away from this particular interaction value thermal- 
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ization occurs on much longer time scales. 

The critical interaction can be identified for example by 
plotting the time evolution of the quantity^ 

An(t) = n(e k = 0_, f) - n(e k = 0+, t), (30) 

which is the size of the discontinuity of the momentum dis- 
tribution function n(ek, t) = (cl a (t)ck,cr(t)) at the Fermi en- 
ergy. After a quench to U < [/^ uench 5 the system initially set- 
tles into a nonthermal quasistationary state characterized by a 
nonzero An(i)^For U > [/ c quench , An(t) exhibits collapse- 
and-revival oscillations. At U w Uc Uench , An(t) rapidly van- 
ishes, with no sign of trapping in an intermediate nonthermal 
state. 

A very similar relaxation dynamics can be observed af- 
ter an interaction ramp with relatively short ramp time r. 
In Fig. [TO] we plot An(t) for different values of the inter- 
action. The top panel shows results for linear ramps with 
ramp time r = 1.25, and the lower panel for ramps of the 
form r(x) = x + 0.87 sin(47ra;) with ramp time r = 1.25. 
The critical interaction strengths (just before the onset of 
the collapse-and-revival oscillations) are 

jjHneai ~ 3 75 and 

^oscillating ^ 4 25. Apparently, the critical interaction strength 
for given ramp time r depends on the ramp shape and this de- 
gree of freedom can be used to design protocols for which the 
desired interaction strength corresponds to a dynamical phase 
transition point. 

The green curves with triangles and the insets of Fig. [TTJ 
demonstrate that also in the ramp case, the critical interaction 
strength is associated with fast thermalization. At the critical 
point, An(t) vanishes within a few inverse hoppings, both for 
linear and oscillating ramps. In the insets we show a compar- 
ison of n(e, t = 3.125) with an equilibrium distribution func- 
tion for the interacting system (temperature T = 0.547 for the 
linear ramp and T = 0.551 for the oscillating ramp). The tem- 
peratures of the thermalized systems have been computed by 
comparing the total energy after the ramp to the temperature- 
dependent total energy obtained from equilibrium DMFT sim- 
ulations. The fast thermalization at the dynamical transition 
and the ramp shape dependence of the critical interaction 
strength imply that it is possible to prepare thermal equilib- 
rium states over a range of interaction values by suitably de- 
signing the ramp protocol. 

V. CONCLUSION 

Using nonequilibrium DMFT, we investigated the excita- 
tion energy of the one-band Hubbard model after ramping up 
the interaction within a time r of the the order of a few in- 
verse hoppings. Based on the perturbative analysis of Ref.l24l 
we determined optimal ramp protocols within a general set of 
piecewise linear ramp functions. This analysis indicates that 
the optimal ramp shape is oscillating around the linear ramp 
U(t)/U(T) = t/r, with a period which is determined by the 
support of the excitation density R(io). (For small values of 
U, the latter depends mainly on the density of states.) Apply- 
ing these optimized ramp shapes to larger interactions (U = 3 
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FIG. 10: Time evolution of the jump An(e, t) in the momentum 
distribution function for indicated values of the ramp amplitude 
(ramp time r = 1.25). The top panel shows results for linear 
ramps, and the bottom panel results for an oscillating ramp r(x) = 
x + 0.87sin(47ra;). For U « 3.75 (linear ramp) and U « 4.25 
(optimized ramp), the jump vanishes within a time of less than 1.5 
inverse hoppings. The insets compare the distribution functions at 
n(e, t = 3.125) to the distribution functions of equilibrium systems 
with interaction U and a total energy identical to the energy after the 
ramp (top panel: U = 3.75, T = 0.547, bottom panel: U = 4.25, 
T = 0.55). 



for a model with bandwidth 4) yields considerable reductions 
in the effective temperature, compared to linear ramps. For 
short ramp times r < 1.5, both unconstrained paths and paths 
constrained to the interaction range [0, U] result in effective 
temperatures which are 30% to 50% lower than those obtained 
with linear ramps. 

Up to intermediate interaction values, the optimal ramps 
determined from the perturbative approach coincide well with 
those determined from the QMC calculations. Our analysis 
thus suggest that one may use the perturbative formula to de- 
termine optimal ramps when an optimization using QMC sim- 
ulations is not possible. This is true in particular for slow 
ramps with ramp times up to 100 inverse hoppings, which 
are used for the preparation of states in cold atom systems. 
Guided by the optimal ramp shape at short times (oscillations 



9 



superimposed on a linear ramp), one can, e.g., try to improve 
on the linear ramp by optimizing the amplitude a in a ramp of 
the form U(t)/U(r) = t/r + asm(uit). Using the "optimal" 
frequency u) from Fig. |9j), for example, the perturbative ap- 
proach predicts that an optimization of a results in a decrease 
of the excitation energy of more than 30% with respect to the 
linear ramp, even in the limit of large r, where the absolute 
value of the excitation energy becomes small (~ 1/r 2 ). 

The precise shape (oscillation frequency) of the optimal 
ramps found in this study is specific to the U = initial 
state. For ramps within the insulating phase, the ramp spec- 
trum will be determined by the interaction energy U rather 
than the bandwidth, and the optimal ramp shapes will be dif- 
ferent. However, in view of our results a perturbative analysis 
should still provide a numerically efficient way to compute 
ramp shapes which result in low excitation energies. 

In addition to minimizing the excitation energy, we have 
briefly addressed the question of thermalization after the 
ramp. We have demonstrated that the thermalization after a 



ramp strongly depends on the final interaction and the ramping 
parameters. Our results indicate the existence of a dynamical 
transition (associated with rapid thermalization), similar to the 
behavior after an interaction quench. The ramp shape depen- 
dence of this dynamical transition point may be exploited to 
design ramp protocols which yield a thermal equilibrium state 
within a time of a few inverse hoppings after the switch-on of 
the interaction. 
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